── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 3.5.0 ✔ tibble 3.2.1
✔ lubridate 1.9.3 ✔ tidyr 1.3.1
✔ purrr 1.0.2
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Código
library(lme4)
Loading required package: Matrix
Attaching package: 'Matrix'
The following objects are masked from 'package:tidyr':
expand, pack, unpack
Código
library(lmerTest)
Attaching package: 'lmerTest'
The following object is masked from 'package:lme4':
lmer
The following object is masked from 'package:stats':
step
Código
library(merDeriv)
Loading required package: nonnest2
This is nonnest2 0.5-6.
nonnest2 has not been tested with all combinations of supported model classes.
Loading required package: sandwich
Loading required package: lavaan
This is lavaan 0.6-17
lavaan is FREE software! Please report any bugs.
Rows: 6945 Columns: 4
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): ISO2
dbl (3): SALID1, YEAR, deaths
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
gbmt_log_merged %>%group_by(Clusters, group, YEAR) %>%summarise(deaths_per_100k =mean(deaths_per_100k) ) %>%ungroup() %>%ggplot(aes(x = YEAR, y = deaths_per_100k, color =as.factor(group))) +geom_line(aes(group = group)) +geom_point() +labs(title ="Promedio de Muertes por 100,000 habitantes por Grupo a lo largo del Tiempo",x ="Año", y ="Muertes por 100,000 habitantes", color ="Grupo") +facet_wrap(vars(Clusters)) +theme_minimal()
`summarise()` has grouped output by 'Clusters', 'group'. You can override using
the `.groups` argument.
Standardizing coefficients only works for fixed effects of the mixed
model.
Model Summary
Parameter
Std. Coef.
SE
95% CI
t(3976)
p
(Intercept)
0.00
0.00
(0.00, 0.00)
-24.26
< .001
YEAR
0.10
4.17e-03
(0.10, 0.11)
25.11
< .001
group (2)
-5.84
1.27
(-8.33, -3.34)
-4.59
< .001
group (3)
10.65
1.26
(8.19, 13.11)
8.48
< .001
group (4)
3.80
1.18
(1.49, 6.10)
3.23
0.001
YEAR * group (2)
5.86
1.27
(3.38, 8.35)
4.62
< .001
YEAR * group (3)
-10.83
1.25
(-13.29, -8.37)
-8.64
< .001
YEAR * group (4)
-3.89
1.17
(-6.19, -1.59)
-3.32
< .001
6.1.2 Poisson
Código
model_3_poisson <-glmer( deaths ~ YEAR * group + (1|SALID1), data = gbmt_log_merged %>%filter(Clusters ==3),family =poisson(link ="log"))
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
- Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
- Rescale variables?
model_4_poisson <-glmer( deaths ~ YEAR * group + (1|SALID1), data = gbmt_log_merged %>%filter(Clusters ==4),family =poisson(link ="log"))
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
- Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
- Rescale variables?
Estimación de los modelos con la información poblacional como offset
Código
model_3_poisson_off <-glmer( deaths ~ YEAR * group + (1|SALID1) +offset(log(population_imp)), data = gbmt_log_merged %>%filter(Clusters ==3),family =poisson(link ="log"))
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 0.00333122 (tol = 0.002, component 1)
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
- Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
- Rescale variables?
model_4_poisson_off <-glmer( deaths ~ YEAR * group + (1|SALID1) +offset(log(population_imp)), data = gbmt_log_merged %>%filter(Clusters ==4),family =poisson(link ="log"))
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
- Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
- Rescale variables?
Debido a que se ha agregado a los interceptos de las ciudades como efecto aleatorio, entonces las estimaciones y predicciones van a variar de ciudad a ciudad. Para poder visualizarlo, se podría resumir los datos con su media.
Código
plot_poisson_offset <- gbmt_log_merged2 %>%group_by(Clusters, group, YEAR) %>%summarise(across(c(predicted_deaths_100k, deaths_per_100k), mean ) ) %>%ungroup() %>%ggplot(aes(x = YEAR, y = predicted_deaths_100k,color = group)) +geom_point(aes(y = deaths_per_100k), alpha =0.5,position =position_jitter(h =0.2)) +# Datos observadosgeom_line() +labs(title ="Promedio de Muertes por 100,000 habitantes por Grupo",x ="Año", y ="Muertes esperadas por 100,000 habitantes", color ="Grupo") +facet_wrap(vars(Clusters)) +theme_minimal()
`summarise()` has grouped output by 'Clusters', 'group'. You can override using
the `.groups` argument.
plot_poisson2_offset <- gbmt_log_merged3 %>%group_by(Clusters, group, YEAR) %>%summarise(across(c(predicted_deaths_100k, deaths_per_100k), mean ) ) %>%ungroup() %>%ggplot(aes(x = YEAR, y = predicted_deaths_100k,color = group)) +geom_point(aes(y = deaths_per_100k), alpha =0.5,position =position_jitter(h =0.2)) +# Datos observadosgeom_line(linewidth =1) +labs(title ="Promedio de Muertes por 100,000 habitantes por Grupo a lo largo del Tiempo",x ="Año", y ="Muertes esperadas por 100,000 habitantes", color ="Grupo") +facet_wrap(vars(Clusters)) +theme_minimal()
`summarise()` has grouped output by 'Clusters', 'group'. You can override using
the `.groups` argument.
plot_poisson2_offset_g <- gbmt_log_merged2 %>%ggplot(aes(x = YEAR, y = predicted_deaths_100k,color = group)) +geom_point(aes(y = deaths_per_100k), alpha =0.5,position =position_jitter(h =0.2)) +# Datos observadosgeom_line(linewidth =1) +labs(title ="Promedio de Muertes por 100,000 habitantes por Grupo a lo largo del Tiempo",x ="Año", y ="Muertes esperadas por 100,000 habitantes", color ="Grupo") +facet_wrap(vars(Clusters)) +theme_minimal()plot_poisson2_offset_g
Debido a que se ha agregado a los interceptos de las ciudades como efecto aleatorio, entonces las estimaciones y predicciones van a variar de ciudad a ciudad. Para poder visualizarlo, se podría resumir los datos con su media.
Código
plot_binomialoffset <- gbmt_log_merged4 %>%group_by(Clusters, group, YEAR) %>%summarise(across(c(predicted_deaths_100k, deaths_per_100k), mean ) ) %>%ungroup() %>%ggplot(aes(x = YEAR, y = predicted_deaths_100k,color = group)) +geom_point(aes(y = deaths_per_100k), alpha =0.5,position =position_jitter(h =0.2)) +# Datos observadosgeom_line() +labs(title ="Promedio de Muertes por 100,000 habitantes por Grupo",x ="Año", y ="Muertes esperadas por 100,000 habitantes", color ="Grupo") +facet_wrap(vars(Clusters)) +theme_minimal()
`summarise()` has grouped output by 'Clusters', 'group'. You can override using
the `.groups` argument.
6.5 Modelo 5 - Con offset sin efectos aleatorios - Binomial Negativo
6.5.1 Estimación
Las estimaciones no consideran variaciones entre las ciudades
Código
model_3_binomial_neg2 <-glm.nb( deaths ~ YEAR * group +offset(log(population_imp)), data = gbmt_log_merged %>%filter(Clusters ==3))performance(model_3_binomial_neg2) %>%print_html()
model_4_binomial_neg2 <-glm.nb( deaths ~ YEAR * group +offset(log(population_imp)), data = gbmt_log_merged %>%filter(Clusters ==4))performance(model_4_binomial_neg2) %>%print_html()
plot_binomial2offset <- gbmt_log_merged5 %>%group_by(Clusters, group, YEAR) %>%summarise(across(c(predicted_deaths_100k, deaths_per_100k), mean ) ) %>%ungroup() %>%ggplot(aes(x = YEAR, y = predicted_deaths_100k,color = group)) +geom_point(aes(y = deaths_per_100k), alpha =0.5,position =position_jitter(h =0.2)) +# Datos observadosgeom_line(linewidth =1) +labs(title ="Promedio de Muertes por 100,000 habitantes por Grupo a lo largo del Tiempo",x ="Año", y ="Muertes esperadas por 100,000 habitantes", color ="Grupo") +facet_wrap(vars(Clusters)) +theme_minimal()
`summarise()` has grouped output by 'Clusters', 'group'. You can override using
the `.groups` argument.
plot_binomial2offset_g <- gbmt_log_merged5 %>%ggplot(aes(x = YEAR, y = predicted_deaths_100k,color = group)) +geom_point(aes(y = deaths_per_100k), alpha =0.5,position =position_jitter(h =0.2)) +# Datos observadosgeom_line(linewidth =1) +labs(title ="Promedio de Muertes por 100,000 habitantes por Grupo a lo largo del Tiempo",x ="Año", y ="Muertes esperadas por 100,000 habitantes", color ="Grupo") +facet_wrap(vars(Clusters)) +theme_minimal()
Rows: 371 Columns: 2
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (2): SALID1, PUBSALID1
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Reading layer `SALURBAL_L1AD_NoSmallIslands_PublicIDs20230922' from data source
`/home/brian/data/Github/salurbal_2/01_data/raw/Data Request MS242_09222023/MS242_L1.gdb'
using driver `OpenFileGDB'
Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : GDAL
Message 1: organizePolygons() received a polygon with more than 100 parts. The
processing may be really slow. You can skip the processing by setting
METHOD=SKIP, or only make it analyze counter-clock wise parts by setting
METHOD=ONLY_CCW if you can assume that the outline of holes is counter-clock
wise defined
Simple feature collection with 371 features and 10 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -117.1243 ymin: -54.28742 xmax: -34.79288 ymax: 32.71865
Geodetic CRS: WGS 84
